Simulation of thermal conductivity and heat transport in solids 
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Using molecular dynamics (MD) with classical interaction potentials we present calculations of 
thermal conductivity and heat transport in crystals and glasses. Inducing shock waves and heat 
pulses into the systems we study the spreading of energy and temperature over the configurations. 
Phonon decay is investigated by exciting single modes in the structures and monitoring the time 
evolution of the amplitude using MD in a microcanonical ensemble. As examples, crystalline and 
amorphous modifications of Selenium and Si02 are considered. 



PACS number(s): 63.20, 63.50, 65.40, 65.70 
I. INTRODUCTION 

Thermal conductivity and heat transport play an impor- 
tant role in the understanding of structural and dynam- 
ical differences between amorphous and crystalline sub- 
stances; however, the underlying mechanisms are not yet 
fully understood. Since the low temperature experiments 
by Zeller and Pohl thermal conductivity and specific heat 
are knowrLto show a universal and anomalous behaviour 
in glassesEI. Phenomenological models to explain these 
outstanding effects have been proposed by P.W. Ander- 
son et alJ3 and by Phillips et aljj (TLS-model: tunnel- 
ing in two- level systems) , and an extension to somewhat 
hi&her temperatures has been worked out by Karpov et 
al.Q in order to take into account contributions caused by 
anharmonic vibrations and thermally activated hopping 
processes or relaxations in glasses (soft potential model) . 
However, a fully consistent picture requires the analysis 
of both equilibrium, steady-state and non-equilibrium as- 
pects of the transport properties of solids. 
During the last two decades non-equilibrium molecular 
dynamics (NEMD) has been applied to study such prop- 
erties, e.g. heat transport, the evolution of shock waves, 
and decay of phononsB tl Based on these pioneering in- 
vestigations it becomes possible to study the properties 
of solids using realistic interaction potentials, enabling us 
to distinguish between the influence of specific features, 
e.g. the nature of chemical bonding, and a more universal 
behaviour, e.g. a more general physical dynamics. Since 
more than -25-jyears shock waves have been simulated us- 
ing NEMDBtj with the focus on shock wave propagation, 
on plastic, deformation of solids or shock wave-induced 
melting O Very recently large-scale computer simulations 
have been performed to gain insight into plastic deforma- 
tions in solids induced by shock waves, where the influ- 
ence of dislocations and defects on-the plastic deforma- 
tion is investigated using NEMD £3 In these extensive 
simulations boundary effects could be excluded. 
The estimation of phonon lifetimes and their influeacej 
on transport properties is stressed by several authors E't^l 
Ladd et al. have calculated the thermal conductivity for 



an FCC lattice using both density- and heat-flux corre- 
lation functions. tl Allen and Feldman have investigated 
the thermal conductivity in silicon by deriving and eval- 
uating a formula-based on the Kubo and Greenwood- 
Kubo formalisrnEj representing the heat current opera- 
tor in terms of eigenmodes determined from the dynamic 
matrix of the disordered Si-structure. 
Using MD-simulations, Michalski has studied the in- 
fluence of harmonic and anharmonic contributions to 
the atomic interaction potential on thermal conductivity 
and heat transport in two-dimensional (quasi-crystalline) 
structures. Furthermore, he has investigated the influ- 
ence of delocalized and localized modes, the latter being 
responsible for strong anharmonic effects in glasses O 
Our work is concerned with molecular dynamics sim- 
ulations of thermodynamical steady-state and non- 
equilibrium transport properties in realistic covalent 
structures, where the main contributions to the ther- 
mal conductivity k and the propagation of heat pulses 
and shock waves originate from the vibrational degrees 
of freedom. i-Qur simulations for k are similar to those 
of MichalskiliJ, but differ from this work by our use 
of three-dimensional structures with periodic boundary 
conditions in all lateral dimensions. 
One aim of our simulations is to check whether one 
can simulate and realistically describe thermodynami- 
cal properties with "usual" interatomic potentials and 
structures in a realistic way, and how reliable such poten- 
tials and configurations are in modelling heat transport 
in solids. Another issue we want to address is the im- 
portance and influence of structural differences in glasses 
and crystals regarding heat transport. Here, the com- 
puter experiment mimicking the macroscopic set-up used 
to calculate k, is complemented by the study of phonon 
decay and wave propagation in such solids. 
In the next section, we briefly describe the systems and 
interaction potentials which we have used. In section 
III, we explain the methods in detail. The results and 
comparison with experiments are presented in section IV , 
followed by a discussion of the results. 
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II. EXAMPLE SYSTEMS 



As representative covalent materials we have used Se and 
Si0 2 . 

Selenium readily ■— ■ forms 

glasses and amorphous structures .113 Several crystalline 
structures exist, including two (a— and (3— ) monoclinic 
forms with four eight-membered rings packed differently 
in the unit cell. Under standard conditions the most sta- 
ble crystalline phase consists of .infinite parallel helical 
chains with trigonal symmetryE_3ii3. We used this so- 
called trigonal (t-)Se for our investigations of crystalline 
selenium. 

The potential used to simulate Se has been described 
earlier .113 The parametrization of the potential was cho- 
sen to mimic certain structural properties of selenium: 
The potential was fitted to reproduce bond lengths—an.- 
gles and bonding energies of small Se-moleculesESEa, 
and to | _aiy | e- 1 a reasonable description of the trigonal 
crystal»El 

SiC>2 exists in many different crystalline allotropes (e.g. 
a- and /3-quartz, high- and lop,cristobalite, tridymite, 
keatite, coesite andp£tisiiovite)E2l and is known to be a 
strong glass former .E3~t3 The atomic interaction poten- 
tial was fitted by Vashishta et al.c3 in order to reproduce 
structural and dynamical properties of both crystalline 
and glassy phases. 

From experiment it is known that a-Si02 has a quite, 
high thermal conductivity compared to others glassesE3 
which have typically thermal conductivities of the order 
of several O.lW/mK. 

Using these empirical potentials for Se and SiC>2 glasses 
are generated, fay . rapid MD-qucnchcs from well equili- 
brated liquidsEiEy and a final quench to K. 



III. METHODS AND CALCULATIONS 



With molecular dynamics (MD) one can_S|imulate and 
investigate properties of complex systems. ETEj In order 
to determine structural, dynamical and thermodynami- 
cal quantities one typically explores correlation functions 
(e.g. Van-Hove-correlation to get an insight into the ra- 
dial distribution of atomic distances, and vel«city-auto- 
correlationEZI or displacement-auto-correlationE3 to calcu- 
late the vibrational spectrum of configurations). Here we 
want to describe a more direct way to determine thermal 
properties of solids. 



A. Thermal conductivity 

Experimentally the thermal conductivity is determined 
by measuring the stationary heat flux necessary to main- 
tain a tempjesature gradient generated by two heaters in 
the sampled One fundamental problem in experiment 
is the thermal insulation of the sample against heat loss 



to the surroundings. Note that in computer experiments 
this problem can be avoided by imposing periodic bound- 
ary conditions. 

In this paper we describe a simulation which is designed 
to mimic the experimental setup. First we equilibrate 
configurations (with N = 1470 to 23520 particles) at tem- 
perature T for several thousand MD time steps (typically 
4000 to 40000 MDS). In order to determine the thermal 
conductivity k of a given structure, we select two '"con- 
tact" ' layers of atoms (typically about 10% of all atoms 
comprising the structure) and fix the average tempera- 
tures in these layers at T L = T± AT by scaling the atomic 
velocities according to the formula ^NksT — Ekm- For 
symmetry reasons these layers are separated by half the 
box-length. In order to reduce surface effects we apply 
periodic boundary conditions in all spatial directions. 
The setup of our computer experiment is sketched in 
Fig. [l]. The atoms outside these layers are allowed to 
move according to Newton's laws. The velocities of these 
particles are not rescaled. After some time (typically 
several thousand MDS) a stationary temperature gradi- 
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has developed. To determine the temperature 

gradient, we calculate the "local" temperature (i.e. the 
kinetic energy) of sub-layers of the structure by divid- 
ing the system into parallel layers of equal thickness in 
which we measure the "local" (kinetic) temperature as 
described in Refs. ||,[40| In Fig. || we plot the tempera- 
ture of the sub-layers versus the z-coordinate of the layers 
(z 6 [—L/2, L/2] corresponds to the original box, the ad- 
ditional points are plotted to show the periodic boundary 
conditions). The fluctuations of the mean temperatures 
in those sub-ensembles for which we do not regulate the 
temperature are about ±4%, whereas in the two layers, 
where the temperature is fixed by scaling the velocity, the 
local temperatures change about ±2% (this is the change 
of the local temperature in one time step before scaling) . 
Note the clear development of a linear temperature gra- 
dient during the simulation. This shows the development 
of a steady non-equilibriurp-S|tate. which enables us to use 
Fourier's law of heat flowllirealill, where the heat flux j 
necessary to maintain the temperature difference 2AT is 
given by: 
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Here, At is the time step used in the MD simulation, A is 
the interface area of the sample perpendicular to the heat 
flux, and (AE) is the average energy per time step At 
which is added and subtracted, respectively, in the layers 
representing the heat contacts. The changes of the energy 
(AE) are determined by the changes of the temperatures 
of the layers: We heat and cool the layers by rescaling 
the velocities of the particles comprising the layers, i.e. 
changing the atomic velocities from Vj (before scaling) 
to (after scaling). Therefore, the energy difference 
is given by (AE) = (AE kin ) = <± m,(vf - v 2 )) = 
|iVifcs (STl), with Nl the numbers of atoms in the layers 
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(large enough to define a sensible local temperature) and 
(STl) the temperature change necessary to maintain the 
desired temperature Tl — T ± AT of the layers. The 
amount (STl) is averaged over the last 20000 MDS of our 
simulation. To avoid quantum effects, which we cannot 
account for (e.g. tunneling in two- level-systems), we have 
chosen "intermediate" temperatures to simulate and to 
measure thermal conductivity. 

In order to test our "computer experiment" , we vary (a) 
the temperature gradients, (b) the layer thicknesses, (c) 
the temperatures, and (d) the system size. 

(a) We find that temperature differences between the con- 
tact layers ranging from 20 to 40% of the average temper- 
ature are sufficient to reach convergence in the "measure- 
ment" of (AE) in reasonable computer time. However, in 
typical computer simulations the temperature gradients 
are of the order 10 10 K/m, which is orders of magnitudes 
higher than in experiments. Such large gradients might 
have a strong influence on the decay of phonons. 

(b) We have varied Nl from 1% to 20% of the total num- 
ber of atoms N in the simulation cell. In our experience 
we get fast convergence of the results by using about 
Nl = 0.1 N atoms in the layers. 

(c) The temperatures used in our simulations range from 
several K upwards to about 400 K. At temperatures be- 
low 1 K, the physical effects in solids are dominated by 
quantum effects which our classical simulation method 
cannot describe. On the other hand, very high temper- 
atures lead to large displacements of the atoms in the 
layers, and the interfaces of the layers and the rest of 
the system roughen considerably. Simulating Se-glass at 
temperatures ranging from 3.5 K to 170 K, we observe 
mean displacements (u) — 0.6 to 0.8 A per atom. In trig- 
onal Se the mean displacement of the particles is about 
0.2 A per atom at T = 370 K. 

The drift in the total energy of the configurations is less 
than 4 parts in 10 4 for the highest temperature employed 
and less than 3 parts in 10 6 for the lowest temperature 
simulated. The reason for this energy drift is the use of 
an isothermal "sub-ensemble" in the part of the system 
where we fix the temperature, while treating the rest of 
the system as a microcanonical ensemble where the total 
energy is conserved. 

(d) From experiment one knows that in "perfectly" crys- 
talline but finite structures the thermal conductivity de- 
pends on the lateral dimensions of the samples. The 
smaller the crystalline sample, the smaller is the mea- 
sured thermal conductivity. In such confined structures 
the mean free path of the phonons is limited by the sur- 
faces of the experimental probes which act as scatterers 
for thc.phonons. This phenomenon is known as Casimir 
limit o In computer simulations where one deals with 
small (typically nano-scale) structures, the influence of 
the system size will be even stronger than in experiments 
with sample diameter of the order of mm. In order to deal 
with this phenomenon, we have performed simulations 
for several system sizes L and extrapolated the thermal 
conductivity for L — » oo. In Fig. we show the inverse of 



the thermal conductivity versus the inverse of the layer 
distance for t-Se. 



B. Heat pulse 

To study the dissipation and spread of kinetic energy in 
the system we induce a heat pulse into the system after 
it has equilibrated. We scale the atomic velocities in a 
layer with Nl particles (typically 15% of the sample) to 
be much higher than the averaged velocities of the atoms 
outside this layer. After this initial heat pulse the sample 
is allowed to evolve freely. Again we calculate the "local" 
temperatures for a sequence of sub-layers. c!3 
The time evolution of the local temperature T(t) in the 
different layers provides insights into the mechanisms un- 
derlying the energy transfer in the systemE In partic- 
ular comparing the evolution of crystalline and glassy 
systems elucidates differences in heat transport caused 
by the structure of the samples. Especially the Fourier 
transform of T(t) can be used to identify modes and vi- 
brations responsible for heat transport. 
We have applied heat pulses to amorphous Se and trig- 
onal Se, both parallel and perpendicular to the helical 
chains. 



C. Shock waves 

The sound velocities of the structures can be determined, 
from the effects of moderate shock waves in the systcm.Eil 
We generate such a wave by shifting the atoms of one 
layer (Nl — 0.02 — 0.05 N) from their equilibrium posi- 
tions by about 2.5% of the nearest neighbour distance in 
one direction and follow the development of the pertur- 
bation. As initial condition the atomic velocities are set 
to zero. Measuring the "local" (kinetic) temperatures [as 
in Ref. ^] in the sub-layers at every time step we derive 
the sound velocity from the space-time evolution of the 
temperature. We have applied this procedure to trigo- 
nal Se both parallel and perpendicular to the chains, to 
a— quartz, and to amorphous Selenium and Si02- 

D. Decay of modes 

One of the central questions in the field of heat transport 
concerns the lifetime of phonons. In perfectly harmonic 
crystals, the lifetime of phonons would be infinite. How- 
ever, due to anharmonic contributions to the potential, 
which gain importance with increasing temperature, the 
mean free path of the phonons is limited by umklapp 
processes. On the other hand, at very low temperatures, 
when the umklapp processes are effectively frozen out, 
the only limiting factors of thermal conductivity and of 
lifetimes/mean free paths of phonons are the scattering of 
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phonons due to surface effects (Casimir limit) and/or im- 
purities or defects of the structures. In real materials the 
thermal conductivity is also limited by electron-phonon 
interactions and isotope effects.c3o 
Clearly, the description of processes in terms of phonons 
depends strongly on whether the exact states in anhar- 
monic systems can be approximated by phonons, i.e., 
whether the phonons decay slowly enough to allow a 
meaningful description of the instantaneous state of the 
system in terms of vibrational modes. 
The phonons, which correspond to the eigenvectors (EV) 
of the dynamical matrix of the system, are determined 
in harrpenic approximation by diagonalization of this 
matrix.ea Since the number of eigenvectors (3N, N be- 
ing the number of particles in the simulation cell) is very 
large, we can only estimate the lifetimes for a subset of 
the EVs. In order to determine the lifetime of a phonon, 
we excite an eigenmode ej (normalized to unity) in the 
system by displacing the atoms according to their con- 
tributions to the corresponding eigenvector, and follow 
the time evolution of the atomic motion using NVE-MD 
where the number of atoms, the volume and the total 
energy are kept constant. The original displacement is 
given by: 



A. Thermal conductivity 



Af(0) 



(2) 



where a is the amplitude of the displacement vector. 
During the MD-simulation we calculate the projection of 
the actual atomic displacement onto a subspace of eigen- 
modes {e m }: 
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\Af(t)\ 
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The projection on the eigenmode Cj will have the pe- 



riod 



1/vj of the excited eigenmode, and the am- 



plitude will be constant as long as no interactions with 
other modes occur, i.e. as long as the mode does not 
decay. Due to the anharmonic interactions between the 
eigenvectors other modes will become excited. The cal- 
culation of the expansion coefficients c^ e (with m ^ j) 
allows to monitor the excitation of other modes e m . With 
this procedure one can follow the time evolution of single 
phonons or modes of the solids. The usual dynamic struc- 
ture factor S(Q, v) can be used to obtain information 
on the broadening of peaks and frequency shifts of typ- 
ically groups of modes, but only little knowledge about 
the types of modes involved and their detailed interac- 
tions. The development of single rnp4es can be extracted 



rie rapt < 
tionD'E. 



in the one-phonon approximationlrEl Using the projec- 
tion procedure, we have investigated the t-Se crystal and 
Se-glasses. 



IV. RESULTS 



Using the algorithms described in Sec. Ill A we have 



simulated the thermal heat conduction K|| (parallel to 
the chains) of trigonal Se. The limitations of the very 
short distance L between the layers kept at fixed temper- 
atures Tl = T ± AT, which is computationally accessi- 
ble, suggest an extrapolation of the distance to infinity. 
The simulation at T = 25 K yields K = 0.033 W/(cmK) 
for N — 23520 atoms. Extrapolating L — > oo gives 
k(L = oo) = 0.072 W/(cm K), which is about a factor 
of 6 lower than the experimental result HI At T = 90 
K even the largest simulated system with N — 23520 
atoms deviates from the experimental value by more than 
a factor 3.8. In the limit of infinite system size, we find 
«(L = oo) = 0.071 W/(cm K) for trigonal Se. This result 
is 30% lower than the experimental result. At T — 370 K 
we extrapolate k(L = oo) = 0.017 W/(cm K) which de- 
viates from the experimental -value (k = 0.0538 W/(cm 
K) measured at T = 400 K.cll At temperatures above 
T = 350 K photons contribute to the thermal conduc- 
tivity which lead to an increase of k.o The discrepancies 
may be due to some short-comings of the potential which 
we use to model Se, especially the low-frequency phonons 
are described insufficiently.E-3 

In the case of Se-glasses, we find a much better agreement 
between the results of our simulations and the experimen- 
tal findings. In Fig. || we plot the thermal conductivity 
versus temperature with double-logarithmic scales. In 
the whole temperature range covered, the deviation hc=. 
tween the theoretical values and the experimental onesca 
never exceeds 20 — 30%. Furthermore, no significant sys- 
tem size effect has been detected in our calculations. 
Our simulations of SiCVglasses show a different be- 
haviour, however. To gain results of thermal conduc- 
tivities comparable with experiment it was necessary to 
construct glasses with L up to 125 A. From these we 
were able to extrapolate the thermal conductivity in the 
limit L — > oo as upper bound of the computed thermal 
conductivity in Si02-glasses. At T = 60 K we find 
K = 0.59 W/(uuK), a result 25% higher than the exper- 
imental value.cla Here one should note that the phonon 
spectrum calculated using the potential for SiOa-given in 
Sec. H overestimates the low-frequency modesEj, while 
the calculated specimm for Se underestimates the low- 
frequency phonons .113 



B. Heat pulse 

In a trigonal Se crystal consisting of 1470 atoms we ex- 
cited in one layer parallel to the z-axis a heat pulse with 
a local temperature which was 5 times higher than the 
temperature in the rest of the system. The layer com- 
prised 210 atoms (14 chains with 15 atoms). We observed 
the spread of the energy in the system. By symmetry the 
direction of the flow was perpendicular to the chains of 



4 



the crystal. The evolution of the "local" temperature 
(i.e. the kinetic energy) of the layer Tjr(i) and of the rest 
of the system T r (t) are shown in Fig. [5] as clotted and 
solid lines, respectively. From Fig. || one can deduce that 
the velocities of the atoms in the layer with the induced 
higher temperature change with a period of about 0.15 
ps (corresponding to a frequency of 6.5 THz) and addi- 
tional lower-frequency modulations in case of crystalline 
trigonal Se. 

The Fourier transforms of Ti >r (t) are shown in Fig. [(]. 
Since T oc v 2 the frequencies of the temperature changes 
are a factor of 2 higher than the corresponding changes of 
the velocities. However, taking this factor 2 into account 
the power spectrum resulting from the Fourier transform 
of T r (t) resesibles in some features the usual density of 
states DOS.E-OH On the other hand, the power spec- 
trum associated with the temperature evolution of the 
excited layer shows strong contributions at low frequen- 
cies (caused by the strong exponential decay of the tem- 
perature of this layer: the "relaxation time" is about 
350 MDS«0.7 ps). These modes reflect the "dissipative" 
character of the evolution of the (kinetic) temperature 
in the layer exposed to the heat pulse. Apart from this 
dominating peak due to the fast decrease in tempera- 
ture at the beginning of the MD-run, we find familiar 
contributions at the high frequency end caused by bond 
stretching modes, and a broad spectrum towards lower 
frequencies .UO 

In another simulation, we excited a heat pulse in a layer 
parallel to the basal plane of t-Se. Again, the energy 
dissipation of the excited layer causes a strong peak to 
appear at the low frequency part of the spectrum. How- 
ever, we again find bond stretching modes clearly devel- 
oped in the power spectrum. The high frequency part of 
the spectrum of excited vibrations in T r (y) resembles the 
one we calculate in the case of perpendicular heat flow, 
but some deviations in the middle of the spectrum occur. 
This part of the spectrum, is assigned to librations and 
bond-bending vibrations .EJ 

We have performed the same simulations for a Se-glass 
consisting of 1470 atoms. We observed a fast (exponen- 
tial) "heat relaxation" of the induced thermal energy, as 
in the case of t-Se. Again, the power spectrum T r {v) ex- 
hibit some similarities with the typical vibrational DOS 
of a-Se. 



C. Shock waves 



scribed in Sec. Ill C , we monitor the development of the 



local temperature of the layers. The sound velocity of 
the system can be derived by plotting the time necessary 
for the perturbation to reach the various layers along the 
system. In Fig. [?] the times are plotted versus the lay- 
ers along the z-axis, yielding a sound velocity v z = 5690 
m/s. This result agrees well with earlier results obtained 
for the longitudinal sound velocity C33 = 5630 m/s.EZI 
After dividing the trigonal Se in layers oriented paral- 
lel to the chains (and perpendicular to the x-axis) and 
displacing the atoms of one layer in x-direction, we again 
measured the speed of the heat pulse, finding an effective 
sound velocity in the x-direction v x = 4600 m/s. Due to 
the symmetry of the trigonal structure, this result is by 
a factor of cos 30° smaller than the longitudinal sound 
velocity vu = 5330 m/s calculated previously O 
To calculate the sound velocity of a-quartz, we inves- 
tigated the time evolution of a perturbation induced in 
the system by displacing the atoms of the central layer 
from their equilibrium positions by about 0.02 A in z- 
direction. The structure comprising N — 2356 atoms was 
divided into 24 layers perpendicular to the z-direction. 
Following the same procedure as before, we calculated 
the sound velocity v z = 8335 m/s. This result agrees 
well with the longitudinal sound velocity C33 = 8250-m/s 
obtained from the elastic constant C33 « 180 GPa.E3 
Next we investigated both a-Se with N = 1470 parti- 
cles and a-Si02 with N = 3456 atoms. We divided the 
structures into 21 and 24 parallel layers, respectively. 
These layers were oriented perpendicular to the z-axis, 
and we induced a perturbation into the glasses by dis- 
placing all the atoms in the middle of the system in pos- 
itive z-direction by about 0.05 and 0.02 A, respectively. 
As before, we observed the time dependence of the in- 
duced shock wave by following the kinetic temperature 
evolution of the various layers. In both cases we did 
not observe any clearcut wave-like perturbation propa- 
gating through the system. Therefore, it was impossible 
to determine reliably, or at least to estimate the sound 
velocities of the glasses from these simulations. This neg- 
ative outcome of the analysis of the behaviour of shock 
waves in glassy and disordered structures contrasts with 
the findings for the crystalline counterparts and might 
be traced back to structural differences, e.g. the lack of 
long-range order in amorphous structures. Thus struc- 
tural defects in glasses might cause the "overdamping" 
of the wave, which is most likely relate d to t he rather 
fast decay of phonons in glasses (cf. Sec. |[VD|) . 



We have determined the sound velocities of solids from 
the analysis of shock waves both for crystals and glasses 
(t-Se, a-quartz, a-Se, and a-Si02). We divided the trig- 
onal Se-crystal with N = 2940 atoms into 42 layers per- 
pendicular to the chains and studied the time evolution 
of a perturbation induced in the system by displacing 
the atoms of one layer of the structure from their equi- 
librium positions by about 0.05 A in z-direction. As de- 



D. Decay of modes 

Exciting single modes in the system allows us to study 
their ti me evo lution. We applied the procedure described 
in Sec. HID to a subset of modes only, since the trig- 
onal Se-crystallite chosen comprised N = 1470 atoms. 
In the low-frequency regime [y < 5THz) we investigated 
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the decay of n(=25) modes. To estimate the influence 
of temperature, we excited the EVs both for the T = 
K-configuration and for a crystalline configuration equi- 
librated at T = 60 K. In Fig. || we plot the projection 
c rCj [cf. Eq. [| versus observation time t. The initially 
excited mode has a frequency v = 1.39 THz and is de- 
localized, its participation ratio being p u — 0.62. The 
envelope of the projection exhibits an exponential decay 
with a time-constant r w 70 ps. The Fourier transform 
gives insight into the frequency shift of the modes due 
to anharmonic effects and into the spectral broadening 
of the mode. The latter effect is also relevant for the 
lifetime of a mode but not so reliable due to the finite 
simulation time. The Fourier transform of the projection 
c rCj shown in Fig. || leads to the spectral density plot- 
ted in Fig. ^| with a maximum at v = 1.36 THz. This 
shift in frequency means that the mode has softened with 
temperatureJia In order to verify, whether the shift in fre- 
quency is due to temperature effects or is an artefact of 
the damping of the mode, we have calculated the Fourier 
transform of cos(W) exp(— t/r) with ui — 1.39 THz and 
withr w 7-10" 10 ,7-10" 11 ,7-10- 12 s. Even for the short- 
est decay time r we find a shift less than 10~ 4 u>. Thus 
we conclude that the observed shift in frequency is really 
caused by the influence of the temperature on the modes, 
in qualitative, agreement with the softening observed in 
experiments 

Due to interactions between the modes additional fre- 
quencies are excited in the system. In order to estimate 
the excitation of additional modes, we have calculated 
the projections of the actual displacement on modes {e m } 
with frequencies similar to the initially excited one. We 
sum the square of the projections on 50 modes. After 
an observation time of about 35 ps, the square of the 
expansion coefficient of the initially excited mode is still 
larger than 0.6, and the sum for the energetically similar 
modes is 0.3. Obviously, the atomic displacements are 
still restricted to a small subspace of EVs. Thus, only 
1% of the eigenmodes of the configuration are needed to 
describe the atomic motions to more than 90%. 
This detailed analysis can be complemented by a calcu- 
lation of S(Q, v) that has the advantage to be directly 
comparable with experiment. However, in general for a 
given Q several phonons contribute to the same S(Q, v). 
Therefore, an unambiguous assignment of the spectrum 
of S(Q 7 v) to single modes is not possible. As an ex- 
ample we show in Fig. [l^ for two values of Q in the 
[0 £ 0]-direction the resulting S(Q : v). The estimate of 
the lifetime by FWHM (full width half maximum) yield 
for the dominant phonon in Q — 27r/a[0 0.2 0] approxi- 
mately 25 ps. But due to the superposition of the peaks it 
is difficult to estimate the lifetimes of the other phonons. 
In particular, there are several phonons which cannot be 
investigated calculating only S(Q, v) because of this Q- 
selectivity. 

In the high frequency regime we followed the time evo- 
lution of a mode just above the gap of the DOS. The 



eigenfrequency of this (localized) mode was = 6.2 THz, 
and it decayed within 2 ps. However, the simulations of 
extended modes in the high frequency end of the DOS 
did not show any appreciable decay of the EVs within 
the observation time. 

For glasses, the behaviour of the modes is totally differ- 
ent. We have diagonalized the dynamical matrix of one 
Se-glass and calculated the EVs. In the low-frequency 
regime there exist beside extended eigenvectors a few 
quasi-localized modes, while in the high-frequency part 
of the spectrum nearly all modes are localized. In the 
regimes of low, medium and high frequencies we have per- 
formed similar simulations as in the case of trigonal Se. 
For a quasi-localized mode (with a frequency v = 0.31 
THz and a participation ratio p v = 0.09) we observe a 
fast decay of the eigenmode. No clear vibration develops 
or can be identified in the glass. During the MD-run, we 
also checked, whether a subset of about 20 modes with 
similar frequencies might be more stable. When analyz- 
ing phonons in the low-frequency regime we found that 
this small subset of eigenmodes (less than 1% of the EV 
of the configuration) described the atomic motions to an 
extent of about 30%. The phonons of the medium and 
high-frequency regime were much less stable and decayed 
within several ps. 

Finally we considered the decay of phonons in the non- 
equilibrium steady state, in order to estimate the influ- 
ence of the strong temperature gradient on the lifetimes 
of the eigenvectors. We excited some phonons in a t-Se 
crystallite with N = 1470 atoms, where two layers were 
kept at fixed temperatures Tl =T± AT, and monitored 
the projection of the actual displacement onto the ad- 
ditionally excited eigenvector. In comparison with the 
lifetime of the same mode excited in a configuration with 
constant temperature, the phonon proved to be much less 
stable: typically about a factor of 5. However, since we 
extrapolate to infinite system size keeping the tempera- 
ture difference between the contacts constant, the effect 
of the strong temperature gradient present during the 
simulations should be largely eliminated when calculat- 
ing k(L = oo). 

V. SUMMARY AND CONCLUSION 

In this paper we have presented simulations of steady- 
state and non-equilibrium thermal heat transport in con- 
densed matter. As examples we used Selenium and Si02, 
both in crystalline and amorphous modifications. 
In the crystalline structures the thermal conductivity 
showed a clear system size dependence (Casimir- limit), 
due to the restrictions of the mean free path of the 
phonons by the lateral dimension of the simulation cell. 
Extrapolation of the system size L to infinity yielded val- 
ues of the thermal conductivity of the same order of mag- 
nitude as the experimental results. Since there are no 
impurities in the system, the phonons are not scattered 
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at structural defects, vice versa this means that the life- 
times of the eigenmodes are limited only by anharmonic 
effects of the atomic interaction potential. 
In glasses the situation is rather different. In the case 
of Selenium glasses, it was sufficient to simulate struc- 
tures comprising 1470 particles in the simulation cell to 
reach convergence of the results for the thermal conduc- 
tivity. This can be understood considering that the amor- 
phous structure, or to be more precise the state of dis- 
order itself, is the reason for the anomalous and univer- 
sal behaviour of thermal heat transport in glasses. The 
phonons decay fast within several ps. But nevertheless, 
the eigenmodes typically interact with energetically sim- 
ilar modest^ which means that the "subspace" of modes 
is relatively stable. Quartz-glass is knownr-tp possess ex- 
tremely large mean free paths of phonons£3; this agrees 
with our observation that it was necessary to construct 
glass configurations with L w 125 A in order to get re- 
sults comparable to the experimental findings. 
It is a striking result of our work that the agreement be- 
tween theory and experiment for thermal conductivity 
of glasses is better than for the corresponding crystalline 
phases. It seems that the influence of the structural disor- 
der is so strong that the element-specific details of the in- 
teratomic potential are less important than in crystalline 
structures. 

Complementing these steady-state computer experi- 
ments, simulations of shock waves and heat pulses in con- 
densed matter were performed. Again one can clearly see 
the influence of structural order or disorder on the results 
and the outcome of the calculations. 
In the crystalline structures it was possible to identify 
waves moving through the system. This is clearly con- 
nected to the relatively high stability of the modes con- 
tained in the induced wave packet. The behaviour of the 
wave packet reflects the elastic properties of the solid, e. 
g. one can calculate the sound velocity, from which one 
can then deduce the corresponding elastic constant. In 
contrast glasses do not show clearly developed propagat- 
ing waves after inducing a heat pulse or a shock wave 
into the structure. Again, the reason for this behaviour 
can be found in the short lifetimes and mean free paths 
of the eigenmodes/phonons of glasses where the disorder 
plays the role of scatterers for the modes. 
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FIG. 3. The inverse thermal conductivity in [W/cmK] -1 
for different system sizes plotted versus inverse system size in 
reduced units at T = 90 K for t-Se. 



FIG. 4. Thermal conductivity K[W/cmK] of Se-glass plot- 
ted versus temperature T in K in a double logarithmic plot. 
o: simulation, +: experimental values from Ref.ca. 



FIG. 5. Evolution of the local temperature of the layer 
Ti(t) (dotted line) and the rest of the crystal T r (t) (solid 
line) in K vs. time in MDS after a hot spot was induced 
into a layer parallel to the chains a trigonal Selenium crystal- 
lite with N = 1470 atoms. The typical period of temperature 
changes is about 77 fs. 



FIG. 6. Fourier transform (power spectrum) Ti(y) (a) and 
f r {v) (b) of T hT {t) of Fig. |. 



FIG. 7. Times in fs needed for the perturbation to reach 
the layers are plotted vs. the distance of the layers in A along 
the direction of propagation. 



FIG. 8. Projection c ro . of the actual displacement onto the 
initially activated eigenvector for a mode with v — 1.36 THz 
in the t-Se vs. time in MDS. The system is equilibrated at 
T= 60 K. The envelopes of the projections are two exponential 
functions. 



FIG. 9. 
Fig. | 



Fourier transform of the projection c rCj shown in 



FIG. 10. Structure factor S{Q,v) in [0£0]-direction for 
two values of Q: solid line Q = 2-7r/a[0 0.2 0] and dashed line 
Q = 27r/a[0 0.3 0]. 



Figure captions 



FIG. 1. Realization of the computer experiment in order to 
determine the thermal conductivity of solids. In all spatial di- 
rections we apply periodic boundary conditions. The shaded 
areas symbolize the layers where the temperatures are fixed. 



FIG. 2. Mean temperatures in the sub-layers of the struc- 
ture vs. z-axis. 
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